arXiv:1501.07884v2 [cs.SY] 10 Feb 2015 


1 


Geometry-based Estimation of Stability Region for 
A Class of Structure Preserving Power Grids 

Thanh Long Vu and Konstantin Turitsyn, Member, IEEE 


Abstract —The increasing development of the electric power 
grid, the largest engineered system ever, to an even more 
complicated and larger system requires a new generation of 
stability assessment methods that are computationally tractable 
and feasible in real-time. In this paper we first extend the 
recently introduced Lyapunov Functions Family (LFF) transient 
stability assessment approach, that has potential to reduce the 
computational cost on large scale power grids, to structure¬ 
preserving power grids. Then, we Introduce a new geometry- 
based method to construct the stability region estimate of power 
systems. Our conceptual demonstration shows that this new 
method can certify stability of a broader set of initial conditions 
compared to the minimization-based LFF method and the energy 
methods (closest UEP and controlling UEP methods). 

I. Introduction 

The electrical power grid is currently undergoing the archi¬ 
tectural revolution with the increasing penetration of renewable 
and distributed energy sources and the presence of millions of 
active endpoints. Intermittent renewable and volatile loads are 
difficult to exactly predict and present challenges concerning 
voltage, frequency, power quality, and power supply during 
unfavorable weather conditions. As such, the existing planning 
and operation computational techniques largely developed 
several decades ago will have to be reassessed and adopted to 
the new physical models in order to ensure secure and stable 
operation of the modern power grids. Among those challenges, 
the extremely large size of the grid calls for the development 
of new generation of stability assessment methods that are 
computationally tractable and feasible in real-time. 

The most straightforward approach to the post-fault sta¬ 
bility assessment problem is based on direct time-domain 
simulations of transient dynamics following the contingencies. 
Rapid advances in computational hardware made it possible to 
perform accurate simulations of large scale systems faster than 
real-time m, El. Alternatively, the direct energy approaches 
El-El, which are accepted and adopted by industry 0, 
allow fast screening of the contingencies while providing 
mathematically rigorous certificates of stability and saving 
more computational resources than time-domain simulations. 
Essentially, the closest UEP method El certifies that the post¬ 
fault dynamics is stable if the system energy at the clearing 
time is smaller than the minimum energy value at every 
unstable equilibrium points (UEP). This method is known 
conservative and not scalable to large-scale power grids since 
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the problem of searching for an exponential number of UEPs is 
an NP-hard problem. The controlling UEP method Q certifies 
that the post-fault dynamics is stable if the system energy at 
the clearing time is smaller than the energy function value at 
the controlling UEP, which is defined as the nearest point on 
the boundary of the actual stability region that the fault-on 
trajectory is approaching, i.e. nearest the fault-cleared state. 
This method is less conservative than the closest UEP method 
since the energy value at the controlling UEP is possibly larger 
than the energy value at the closest UEP. However, as the 
actual stability region is unknown, the controlling UEP can 
only be searched by some heuristic algorithms. 

Recently, we introduced the Lyapunov Functions Family 
(LFF) approach to alleviate some of these drawbacks El. 
The principle of this approach is to provide transient stability 
certificates by constructing a family of Lyapunov functions, 
which are generalizations of the classical energy function, 
and then find the best suited function in the family for given 
initial states. Basically, this method certifies that the post-fault 
dynamics is stable if the fault-cleared state stays within a 
polytope surrounding the equilibrium point and the Lyapunov 
function at the fault-cleared state is smaller than the minimum 
value of Lyapunov function over the flow-out boundary of that 
polytope. Generally, the LEE approach can certify stability of a 
broader set of initial conditions compared to the closest UEP 
method. Also, the introduced optimization-based techniques 
for constructing stability certificates are scalable to large- 
scale power grids, since they avoid identifying the exponential 
number of UEPs. In addition, the LEE approach is applicable 
to stability assessment of power grids with losses 0, which 
is impossible by the standard energy method. 

In this paper, we improve the LEE transient stability assess¬ 
ment method and make two contributions. The first contribu¬ 
tion is the extension of LEE method to structure-preserving 
power systems. The second contribution is a new geometry- 
based method to construct the estimate of stability region of 
the desired equilibrium point, which we argue to possibly be 
larger than that defined by the existing methods. We observe 
that among all of the UEPs, there are many points that are far 
from the equilibrium point and thus are not necessary to be 
counted when we search for the controlling UEP. Therefore, 
we define 2\8 \ points that are the minimum points of Lyapunov 
function over the 2\£\ flow-out boundary segments of the 
considered polytope. Here, \8\ is the number of lines in the 
grids. These 2\8\ minimum points play the role of all possible 
controlling UEPs of the system. The post-fault dynamics 
is certified stable if the fault-cleared state stays within the 
polytope and the Lyapunov function at the fault-cleared state 
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is smaller than the Lyapunov function at the (controlling) 
minimum point corresponding to the polytope’s subset con¬ 
taining the fault-cleared state. This method is less conservative 
than the original LFF method in ||8l since the Lyapunov 
function at the controlling minimum point is possibly larger 
than the minimum value of Lyapunov function over the flow- 
out boundary. In comparison to the controlling UEP method 
we note that since the 2\£\ minimum points play the role 
of all possible controlling UEPs of the system, the proposed 
geometry-based method can certify stability for points for 
which the controlling UEP method cannot. Eurthermore, the 
construction of the minimum points is mathematically rigorous 
and does not involve any heuristic algorithm. Also, knowledge 
of the fault-on trajectory is not required as in the controlling- 
UEP method Q. 

We note that there are many works on Lyapunov function- 
based stability of structure preserving power systems Eol- 
(na. However, the Lyapunov function in these works is usually 
used to prove the local stability of the system; it is not fully 
exploited to construct the stability region of the system as 
in this paper. Instead, in these works the stability region is 
estimated by the energy method. 

11. Structure Preserving Power Systems 

In normal conditions, power grids operate at a stable equi¬ 
librium point. Under some fault or contingency scenarios, 
the system moves away from the pre-fault equilibrium point 
to some post-fault conditions. After the fault is cleared, the 
system experiences the transient dynamics. This work focuses 
on the transient post-fault dynamics of the power grids, and 
aims to develop computationally tractable certificates of tran¬ 
sient stability of the system, i.e. guaranteeing that the system 
will converge to the post-fault equilibrium. In this paper, we 
address this question on a traditional swing equation dynamic 
model of power systems, which is named structure-preserving 
model originally introduced in mni. This model naturally 
incorporates the dynamics of rotor angle as well as response 
of load power output to frequency deviation. However it does 
not model the dynamics of voltage in the system which is the 
main downside of the approach. However, in comparison to 
the classical swing equation with constant impedance loads, 
the structure of power grids is preserved in this approach. 

Assume that the grid has m generators and no buses in 
which no — m buses have loads and no generation. It is 
convenient to introduce fictitious buses representing the in¬ 
ternal generation voltages. So, in the augmented grid we have 
n = no + m buses. Assume that the grid is lossless. The m 
generators have perfect voltage control and are characterized 
each by the rotor angle Sk and its angular velocity Sk- The 
dynamics of generators are described by a set of the so-called 
swing equations: 

mkSk -f dkSk + Pek - Prak = 0, fc = 1,.., m, (1) 

where, mfe is the dimensionless moment of inertia of the gener¬ 
ator, dk is the term representing primary frequency controller 
action on the governor. Pm^. is the effective dimensionless 
mechanical torque acting on the rotor and Pe^. is the effective 
dimensionless electrical power output of the generator. 


Let Pdf. be the real power drawn by the load at bus, 
k = m + 1. .. ,n. In general Pd,. is a nonlinear function 
of voltage and frequency. Eor constant voltages and small 
frequency variations around the operating point P^^, it is 
reasonable to assume that 

Pdk = Pdk + dkdk, k = m + lf...,n, (2) 

where dfc > 0 is the constant frequency coefficient of load. 
When dk = Owe have the constant load model. The electrical 
power Pgfc from the bus into network, where k = 1,..., n, 
is given by 

Pe. = E sin(4 - 4)- (3) 

j^M'k 

Here, the value 14 represents the voltage magnitude of the k*^ 
bus which is assumed to be constant. Bkj are the (normalized) 
susceptance between k*^ bus and bus. Afk is the set of 
neighboring buses of the k*^ bus. Let akj = VkVjBkj- Einally, 
the structure-preserving model of power systems is obtained 
as: 


T dkSk T ^ ^ dkj sin(b/i; Sj) — Pm^f 

k = 1,..., m, 

dkh + E sin(4 - Sj) = - P^^, 

k = m + 1,.. 


(4) 

(5) 


The system described by equations (|4]i-(|5ll has many station¬ 
ary points with at least one stable corresponding to the desired 
operating point. Mathematically, this point, characterized by 
the rotor angles S* = ..., <5*, 0,..., 0]^, is not unique since 

any shift in the rotor angles [<5* -fc,..., b* +c, 0,..., 0]^ is also 
an equilibrium. However, it is unambiguously characterized by 
the angle differences — S* that solve the following 

system of power-flow like equations: 

E o,kj sin{5lj) = Pk,k = l,...,n, (6) 

j&Afk 

where Pk = P^^ ,k = 1,... ,ni, and Pk = -P^^ ,k = m + 
l,...,n. Then, the set of swing equations (|4li-(|5ll is equivalent 
with 


mkdk + dkdk = - E akj{sm{6kj) - sinjb^)), (7) 

jeA'fc 

/c = 1,..., m, 

dkSk = - E afci(sin(4j) - sm(b4)), (8) 

jeA'fc 

k = m + 1,... ,n. 


Eormally, we consider the following problem. 

Transient stability assessment problem: Determine 
if the post-fault scenario defined by initial conditions 
{4(0), 4(0)})!^]^ of the system (|7]l-® leads to the stable 
equilibrium point S* = ..., <5*, 0,..., 0]^. 

We will address this problem by estimating the stability 
region of the stable equilibrium point S*, i.e. the set of 
points from which the system (|7]i-® will converge to S*. 
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Fig. 1. Bounding of nonlinear sinusoidal coupling (sin^^j — sinS^^) by 
two linear functions of angular difference Skj as described in {nj 


If the initial state xq belongs to this estimate set, then the 
corresponding post-fault scenario is determined stable. We 
will use a sequence of techniques originating from nonlinear 
control theory that are most naturally applied in the state space 
representation of the system. Hence, we view the multimachine 
power system (111)-® as a system with the state space vector 
X = [xi,X 2 , x^Y' composed of the vector of generator’s angle 
deviations from equilibrium xi = [Si — S*Sm — ^ 

their angular velocities X 2 = [Si, ... , 6 mY, and vector of 
load’s angle deviation from equilibrium X 3 = [Sm+i — 
S^j^i ,..., —(5*]^. Let E be the incidence matrix of the cor¬ 

responding graph, so that E[5i... SnY = [Yk-Sj){k,j}GeY■ 
Consider matrix C such that Cx = E\5i ... Consider the 
nonlinear transformation F in this representation is a simple 
trigonometric function ^((70::) = [(sin^^j — 

In state space representation the system can be expressed in 
the following compact form: 


Xl = X2 

±2 = M-\-niX 2 - SiE^SF(Cx)) (9) 

±3 = -D2^S2E^SF(Cx) 

where S = 5*1 = [Irnxm ^mxn—mli *^2 ~ 

[^n—mxn—m mxm]- Equivalently, 

x = Ax-BF{Cx), (10) 


with the matrices A, B given by the following expression; 


A = 


Omxm 

Ojnxm 

OmXm 


^mxm 

Om 


Di 


Omxn—m 

Omxn—m 

^mxn—m 


( 11 ) 


and B = [ 0 ^^\s\ -M:[^SiE'^S -DYS 2 E'^S ]^. 

Here, \S\ is the number of edges in the graph defined by the 
susceptance matrix, or equivalently the number of non-zero 
non-diagonal entries in Bkj. 


HI. Lyapunov Functions Family Approach 

This paper proposes a family of Lyapunov functions to 
certify the transient stability for the structure preserving power 
system (fTOl i. The construction of this Lyapunov functions 
family is based on the linear bounds of the nonlinear couplings 


which are clearly separated in the state space representation 
(doll. From Fig. [T] we observe that 


0 < (4, - 4,)(sin4, - sin4^.) < (4, - (12) 


for any |4j +4jl — Therefore, the nonlinearity F{Cx) can 
be bounded by the linear functions in the polytope V defined 
by the set of inequalities |4j + l — ’’’■ 

Exploiting this nonlinearity bounding, we propose to use the 
convex cone of Lyapunov functions defined by the following 
system of Linear Matrix Inequalities for positive, diagonal 
matrices K,H of size 2\8\ x 2\£\ and symmetric, positive 
matrix Q of size 2n x 2n : 


A^Q + QA R 
-2H 


(13) 


where R = QB — C"'"H—{KCAY- For every pair Q, K satis¬ 
fying these inequalities the corresponding Lyapunov function 
is given by 

-'^K{k,j}{cosSkj +dkjSm5lY- ( 14 ) 

Here, the summation goes over all elements of pair set 
8 , and denotes the diagonal element of matrix K 

corresponding to the pair {fc, j}. 

Similar to Appendix A in ®, we obtain the derivative of 
Lyapunov function V(x) along (fTOl i as: 

V(x) = -0.5(Xx - YFY{Xx - YE) - {Cx - FYHE 

= -0.5(Xx - YFY{Xx -YF)-J2 H{k,j}g{k,j}, (15) 


where g{k,j} = {Skj - Sl^ - (sin4i - sinb^j)) (sin4j - 
sinSlj)- From Fig.[T] we have g{k,j} > 0 for any |4j + 4iI ^ 
TT. Hence, V(x) < 0, Vx € V, and thus the Lyapunov function 
is decaying in V. Therefore, we have the following result. 

Theorem 1: In the polytope V, the Lyapunov function de¬ 
fined by Gil is decaying along the trajectory of Cl, i.e., 
V{x{t)) is decaying whenever x{t) evolves inside V. 


IV. Geometry-based Stability Certification and 
Contingency Screening 

A. Construction of Stability Certificate 

In ®, the stability certificate is constructed by finding the 
minimum value I4iin of the function V {x) over the union of 
flow-out boundary segments of the polytope V. Accordingly, 
if the Lyapunov function at the initial state, which stays 
inside V, is smaller than V„iini then the system trajectory is 
guaranteed to converge from the initial state to the desired 
stable equilibrium point. In this paper, we will introduce a 
geometry-based approach for stability certificate construction, 
in which we inscribe inside the polytope V an invariant set 
TZ which is the largest set formed by combining the flow- 
in boundary of the polytope V together with the patches of 
Lyapunov function’s sublevel sets that are guaranteed do not 
meet the flow-out boundary of V. 

In deed, we divide the boundary dVkj of V corresponding 
to the equality [Skj + 4jl = subsets dV^ and 

dV^Y' The flow-in boundary segment dVjfj is deflned by 
[Skj + 4^1 = ^ Skj Skj < 0, while the flow-out boundary 
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Fig. 2. Comparison between the stability region estimates defined by V^min- 
based method and geometry-based method with the stability region obtained 
by the closest UEP energy method (black solid line). The stability region 
estimated by the Vmin method is the intersection of the Lyapunov level set 
(blue solid line) and the polytope defined by —tt — (5* < 5 < tt — (5*. The 
stability region estimated by the geometry-based method is the inner of the 
set whose boundary is combined of solid blue segments 

segment i9'P£“‘ is defined by \5kj + Slj \ = tt and SkjSkj > 0. 
Since the derivative of 6 ^^ at every points on dV^ is negative, 
the system trajectory can only go inside V once it meets 
We define the following minimum values of V(x) on the flow- 
out boundary segment 

where is the flow-out boundary segment of polytope 

V that is defined by 5kj + = ±7r and SkjSkj > 0. Let 

be the point on such that 

Now we consider the set TZ formed by the combination of 
the flow-in boundary of the polytope V together with 

2\£\ segments of Lyapunov function’s sublevel sets. Each of 
these segments goes through one of the 2 \E\ points x°^^^ 
and lies in the half of the polytope V corresponding to 
sign((jfcj) = ±. The conceptual demonstration of the set TZ 
is given as the combination of solid blue lines in Fig. |2l Note, 
these segments can only meet the boundary of T’ at the point 
with SkjSkj = (ttt — < 0, i.e. the point on the flow- 

in boundary. Therefore, the boundary of TZ is composed of 
segments which are parts of Lyapunov function’s sublevel sets 
or flow-in boundary. 

From the decrease of Lyapunov function inside 'P (Theorem 
[T]) we note that from any initial state inside TZ the system 
trajectory cannot escape TZ through the Lyapunov function’s 
sublevel sets. Also, once the system trajectory meets the flow- 
in boundary, it can only go inside the polytope 'P. So, if the 
set TZ is closed, then its inner is an invariant set. In Appendix 

IVlll-Al we prove the following main result of this paper. 

Theorem 2: If the set TZ is c/oiec/Q then the inner ofTZ is 
an estimate of the stability region of the equilibrium point S *, 
i.e., from any initial state xq in the set TZ, the system trajectory 
xt of (nni) will converge to S*. 

* We conjecture that there are always some Lyapunov functions in the family 
defined by the LMIs 1131 such that the set 7Z is closed. In the conceptual 
demonstration of 2-bus system, it is easy to search for such Lyapunov function 
by the adaptation algorithm introduced in O- 


Theorem |2] provides a geometry-based estimate of the 
stability region of the stable equilibrium point. As a conceptual 
illustration, we can observe from Fig. |2] that in the most simple 
case of 2-bus system, the geometry-based method results in the 
largest stability region estimate compared to the closest UEP 
method and the Umin method in |0. 

B. Direct Method for Contingency Screening 

In this section, we will apply the geometry-based stability 
certificate to the contingency screening problem. Essentially, 
the post-fault dynamics is certified stable if the fault-cleared 
state Xo stays within the polytope P and the Lyapunov function 
at Xo is smaller than the Lyapunov function at the (control¬ 
ling) minimum point corresponding to the polytope’s subset 
containing the fault-cleared state. Indeed, for a given fault- 
cleared state Xo, which is determined by integration or other 
techniques, the value of V{xo) can be computed by direct 
application of (fT4l) . If xo is inside the polytope V, we calculate 
the frequency differences Skj. From the \S\ signatures of these 
frequency differences, we can determine the subset of the 
polytope V in which every points have the same signatures for 
frequency differences with Xo ■ Then, from the formulation (fTbl) 
we can define \£\ minimum values V^- , in which is 

either or according to the signature of Skj ■ The 

value of Lyapunov function at the initial state xo should be 
then compared to the minimum of these \£\ minimum values 
If Vo is smaller than this minimum value, the post¬ 
fault dynamics is certified stable, because xo belongs to the 
stability region estimate T* *. 

We note that unlike energy based approaches, the LFF 
method provides a whole cone of Lyapunov functions to 
choose from. This freedom can be exploited to choose the 
Lyapunov function that is best suited for a given initial 
condition or their family. Essentially, we can apply the similar 
iterative algorithm in 0 (Section IV) to identify the Lyapunov 
function that certifies the stability of a given initial condition 
Xo whenever such a Lyapunov function exits. 

V. Simulation Results 

To illustrate the effectiveness of the LFF and geometry- 
based approach in estimating the stability region of power sys¬ 
tems, we consider the classical 2-bus with easily visualizable 
state-space regions. This system is described by a single 2-nd 
order differential equation 

-I- c?i5 -I- asinb — P = 0. (17) 

For this system 5* = arcsin(P/a) is the only stable equi¬ 
librium point (SEP). For numerical simulations, we choose 
m = 1 p.u., d = 1 p.u., a = 0.8 p.u., P = 0.4 p.u., and 
i5* = 7r/6. Figure |2] illustrates the construction of stability 
region estimate for the most simple 2-bus system by the closest 
UEP method, the Umin method in 0, and the geometry-based 
method. It can be seen that there are many contingency scenar¬ 
ios defined by the configuration xo whose stability cannot be 
certified by the energy method, but can be ensured by the LFF 
method. Also, the geometry-based method provides a better 
stability region estimate compared to the Umin method. 
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We can also see that the two minimum points are all 

the UEPs of the system. Hence, the estimate set TZ covers 
the Lyapunov function’s sublevel sets that go through the 
UEPs. Therefore, the geometry-based stability certificate can 
assess transient stability for every initial states in V that the 
controlling UEP method in Q does. 

VI. Conclusions and Path Eorwards 

This paper extended the recently introduced LEE approach 
to transient stability certification of structure-preserving power 
systems. A new geometry-based technique was also introduced 
to further enlarge the estimate of stability region compared to 
the original LEE method. The new estimate is the largest set 
formed by combining the flow-in boundary of the polytope in 
which the Lyapunov function is decreasing together with the 
patches of sublevel sets that are guaranteed do not meet the 
flow-out boundary of that poly tope. Our numerical simulations 
showed that this new estimate of stability region is broader 
than that obtained by the energy methods and the original 
LEE method. In the applications to contingency screening, 
the geometry-based technique in this paper resulted in a more 
complicated algorithm compared to the original LEE method 
in However, the larger stability region estimate obtained by 
the geometry-based method guaranteed that more contingency 
scenarios are screened and certified stable. 

Toward the practical applications of the Lyapunov Eunctions 
Eamily approach to transient stability certification, further 
extensions should be made in the future where more com¬ 
plicated structure-preserving models of power systems are 
considered, e.g. the dynamics of generators’ voltage or effects 
of buses’ reactive power is incorporated in the model. Since 
the LEE method is applicable to lossy power grid m, it is 
straightforward to extend the method to incorporating reactive 
power, which will introduce the cosine term in the model 
Q. This can be done by extending the state vector x and 
combining the technique in this paper with the LEE transient 
stability techniques in ||9l for lossy power grids (without 
reactive power considered). Also, we can see from the proof of 
Theorem [T] that, in order to make sure the Lyapunov function 
is decreasing in the polytope V, it is not necessary to restrict 
the nonlinear terms F{Cx) to be univariate. As such, we can 
extend the LEE method to power systems with generators’ 
voltage dynamics in which the voltage variable is incorporated 
in a multivariable nonlinear function F. 

We envision to develop a new security assessment toolbox 
for practical power grids based on the LEE approach. This tool 
can certify transient stability for a broad set of contingency 
scenarios when the dynamics of power systems in described 
by a number of models, from simple classical reduction model 
to complex structure-preserving model with dynamic voltage 
and reactive power incorporated. Also, this security assessment 
toolbox can certify stability for rather complicated situations 
when the system parameters are changing or unknown via the 
robust stability certificate developed in El. We will build a 
library of models and contingency scenarios the stability of 
which can be certified by this security assessment toolbox. 
This will help us quickly assess the transient stability of 
dynamical power systems by offline algorithms. 
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VIII. Appendix 

A. Proof of Theorem |2]/or Stability Region Estimate 

Since inner of TZ is an invariant set we have x{t) € TZ CP 
for all f > 0. By Theorem [T] we have V{x(t)) < 0 for all t. 
Erom LaSalle theorem, we conclude that the system trajectory 
x{t) will converge to the set {x :V{x)= 0}. This together 
with (ESll imply that the system trajectory will converge to 
the stable equilibrium point S* or to some point lying on the 
boundary of P. However, by the construction of P the second 
case cannot happen. Therefore, the system will converge to 6 *. 
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